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ABSTRACT 

Current X-ray missions are providing high-quality X-ray spectra from neutron stars (NSs) in quiescent low- 
mass X-ray binaries (qLMXBs). This has motivated us to calculate new hydrogen-atmosphere models, in- 
cluding opacity due to free-free absorption and Thomson scattering, thermal electron conduction, and self- 
irradiation by photons from the compact object. We have constructed a self-consistent grid of neutron star 
models covering a wide range of surface gravities as well as effective temperatures, which we make available 
to the scientific community. 

We present multi-epoch Chandra X-ray observations of the qLMXB X7 in the globular cluster 47 Tuc, which 
is remarkably nonvariable on timescales from minutes to years. Its high-quality X-ray spectrum is adequately 
fit by our hydrogen-atmosphere model without any hard power-law component or narrow spectral features. If 
a mass of 1.4 M Q is assumed, our spectral fits require that its radius be in the range R as = 14.51} ■* km (90% 
confidence), larger than expected from currently preferred models of NS interiors. If its radius is assumed to 
be 10 km, then a mass of M ns = 2.2(£g jg M Q is required. Using models with the appropriate surface gravity for 
each value of the mass and radius becomes important for interpretation of the highest quality data. 
Subject headings: radiative transfer — binaries : X-rays — globular clusters: individual (NGC 104) — stars: 
neutron 



1. INTRODUCTION 

One of the primary goals of neutron star (NS) studies is to 
constrain the behavior of matter at high densiti es by measur- 
ing NS masses and radii ( Lattimer & Prakash 2001). Mass 
measurements of high accuracy for a number of radio pulsars 
in close binary systems (often with inferred NS companions) 
are consi stent with a range of NS mass es between 1.25 and 
1.45 M^ dThorsett & Chakrabartvi ri999). Fundamental con- 
straints on NS interior structure can be achieved by mea- 
surement of the gra vitational redshift from the NS surface 
( Cotta met al.ll20 02). Finally, it should be possible to derive 
constraints on the radius of NSs from spectral fits to their X- 
ray emission if the temperature, composition of atmosphere, 
and distance to a NS are known, and the magnetic field is 
sufficiently weak so as not to affect the opacity, or tempera- 
ture distribution, on the NS surface. These requirements can 
be fulfilled for X-ray observations of quiescent low-mass X- 
ray binaries (qLMXBs) containing NSs, particularly those lo- 
cated in globular clusters where the dis tance is well-known 
jBrown et al.lll998l Ikutledge et al. 2002a). In this paper, we 
perform the most accurate such test currently possible, us- 
ing accurate hydrogen atmosphere models constructed specif- 
ically for this project, and a long Chandra observation of a 
globular cluster containing a relatively bright and remarkably 
constant qLMXB. 

Several low-mass X-ray binaries which have been identified 
during outbursts as ac creting NS systems hav e been observed 
in qui escence (see ICampana et alJ 119981 iRutledge et al.l 
2002b). Their quiescent appearance generally differs from 
that of quiescent syst ems cont aining black holes both in their 
X-ray luminosity (Garcia et al.ll200 1[) and their observed X - 
ray spectrum IRutledge et al.Hl 999t iMcClintock et al JI2004D. 
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both of which indicate the presence of a compact object sur- 
face for NS systems and not for black hole systems. These NS 
systems generally show soft spectra, consisting of a thermal, 
blackbody-like component, and possibly a harder component 
extending to higher energies, usually fit with a power-law of 
photon index 1-2 (alth ough some systems are dominated b 
the harder component; Campana et al. 2002; Wiinands et 
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2005a). The thermal component, if fit by a blackbody, pro 
duces inferred radii too small for theoretical NS size esti- 
mates. However, an accreting NS will develop a pure hy- 
drogen atmosph ere if the ac cretion rate falls below ~ 10~ 13 
M yr _1 (iBrown et alJ 119981) . because the metals se t tle ou t 
of the atmosphere within a fe w seconds (Romani 119871) . 
Raiag opal & Romanil (119961) and lZavlin et all (119961) showed 
that hydrogen NS atmospheres shift the peak of the emitted 
radiation to higher frequencies due to the strong frequency 
dependence of free-free absorption. 

A major unsolved question about qLMXBs is the na- 
ture of the X-ray emission. IBrown et alJ (1 19981) advanced 
the idea (also discussed by Campana et alJ Il998h that the 
soft thermal component seen in these field systems can 
be explained by the release, over long timescales, of heat 
injected into the deep crust by pycnonuclear reactions 
driven during accretion (the "deep crustal heating" model). 
This scenario generally predicts the quiescent thermal lu- 
minosity of many qLMXBs, based on their outburst his- 
tory, reasonably well ( Rutledge et al. 2001 a), but not for 
all qLMXBs (cf. IColpi et alJ 120011 bampana et all l2002t 
Wiina nds et alJl2005al) . Thedeep crustal heating model can- 
not explain the hard power-law component, which is often 
attributed to continu ed accretion and/or a shock from a pul- 
sar wind dCampana et alJH998l iBogdanov et al J 120051) . The 
deep crustal heating model also cannot explain the short- 
timescale (~ 10 4 s) v ariability obse rved from Aquila X- l 
JRutledge et alJl2002bl) and Cen X-4 (Cam pana et alJ l2004). 
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Continued accretion has been suggested as an explanation for 
the thermal component, as the radiation spectrum from matter 
accreting radially onto a neutron star should be similar to that 
expected from deep crustal heating (Zampieri et a D 19951) . If 
an absorption feature due to metals in the NS atmosphere 
were to be confirmed in a qLMXB spectrum (as suggested in 
iRutledee et alJl2002bl) . this would provide evidence for con- 
tinued accretion at rates sufficient to explain most or all of the 
thermal emission. 

Applying the hydrogen-atmosphere NS models of 
I Zavrin et alJ (119961) to Chandr a observations of qLMXBs, 
Rutledge et al. (1999, 2001a b) have shown that the radius 
predictions of the models are consistent with the range of 
radii expected from NSs. These analyses have suffered 
from uncertainties in the distances to qLMXB systems, 
and from uncertainties due to fitting two components (the 
power-law plus thermal components) to a spectrum. The 
distance is tightly constrained for so me globular c lusters 
(4% distance uncertainty to 47 Tuc, Grat ton et alJ 12003). 
making the m an excellent target f o r such studies , as pio - 
neered by iRutledge et alJ <2002al) : iGendre et all (|2003b); 
iHeinke et alJ ( I2003IK The well-studied globular cluster 47 
Tucanae (NGC 104; hereafter 47 Tuc) is especially ideal 
due to its close distance (4.85±0.18 kpc), lo w reddening 
(E(B-V) = 0.024 ±0.004, Gratton et alJl2003l) presence of 
two r easonably bright qLMXBs (X5 and X7; Heinke et al. 
120031 her eafter HGL03), an d deep Chandra observations 
(300 ksec. lHeinke et alJ20 05). 

To match the quality of the best Chandra data we desire 
highly accurate hydrogen-atmosphere NS models. We have 
been troubled by the disagreement between the predictions 
of the currently avai lable models of Zavlin et al. ( 199a) an d 
iGansicke et alJ ll2002l) . We also wished to verify whether 
the variation in surface gravity over the relevant range in NS 
mass and radius has a significant effect on the atmosphere 
models and spectral fitting. For these reasons we have pro- 
duced grids of new hydrogen-atmosphere models. We con- 
sider only models of pure hydrogen and for which the mag- 
netic field is sufficiently weak (B < 10 s G) that it may be ig- 
nored in determining the spectrum. The former assumption 
is consistent with previous observations of qLMXBs, as iron 
or solar-abundance atmospheres would be easily identifiable 
by their different spectral shapes (however, subtle departures 
from pure H might still go unnoticed, and affect our results). 
The latter assumption is consistent with the lack of observ- 
able milliseco nd time variability in qLM XBs similar to X7 
(e.g. Aql X-l, Chan dler & Rufledgel2000l) : the accreting mil- 
lisecond X-ray pulsa rs appear to be substan tially fainter and 
harder in quiescence ( Wiinands et al. 2005b). 

A number of codes exist to compute the spectrum for the 
simple case of a pure hydrogen NS atmosphere with no mag- 
netic fi£U i _Sojne_offliese are hmited to zero magnetic field, 
e.g.JRaiagopal & Romanll i ll996l) . lGansicke et al.l(l2002l) . and 
McClinto ck et alJ (12004). Others include magnetic fields, but 
can be applied in the zero field limit, e.g.. lZavlin et alJ (1996) 
and lLlovdl (120031) . 

Even with these simplifying assumptions, grids of NS 
atmosphere models over a wide range of parameters 
are not widely available for use in XSPEC. Some re- 
searchers have computed their own model atmospheres to 
be used in XSPEC and have included rang es of grav- 
ity as well as effective temperature, e.g ., IZavlin et al.l 
(1998), using NSA, and IStage etafl (120041) . using ATM. 
However, the available models included in the standard 



XSPEC package, na mely NSACg/) (jZavlin et alJl!996l) and 
HYD_SPECTRA(g. s ) (Gansic keet alJl2002l) cover a range of 
effective temperatures, but only for the single surface gravity 
logg., = 14.385, corresponding to a "standard" NS model with 
mass M ns = 1.4M© and radius R ns = 10 km [we shall empha- 
size this using the qualifier "(#*)"]• However, as the results 
of this paper will demonstrate, the surface gravity can play an 
important role in the spectral fitting and needs to be taken into 
account. 

A large part of the present effort was devoted to the con- 
struction of grids of models to fulfill this need. The need for 
NS models covering an extensive set of parameters, includ- 
ing surface gravity, motivated us to adapt the code N SAT- 
MOS. This code was previously used by iMcClintock et alJ 
(2004) to investigate the hydrogen atmospheres of hypothet- 
ical, compact objects, some so compact that they lay within 
their own photon spheres; this required taking account of the 
self-irradiation of the surface due to gravitational bending of 
rays. The same model assumptions described in the Appendix 
of McClintock et al. (2004) apply here, with some improve- 
ments: the code now includes the ra diation force in the equa- 
tion of hydrostatic equilibrium (e.g. iPavlov et all 199 II) . and 
includ es heat conduction by elect rons in the energy equation 
(e.g. [kaiagopal & Romanli ri996l) . These and other technical 
improvements made to NS ATMOS are discussed in Appendix 
A. 

As it turned out for the present problem, neither electron 
conduction nor radiation force playe d any substantial role , 
in agreement with the results of, e.g.,|Qansicke et al. (2002). 
Also, the range of parameters where self-irradiation occurs 
does not overlap with the range allowing neutron stars which 
obey causality, indicating self-irradiation is unlikely to be rel- 
evant for neutron stars. On the other hand, allowing for varia- 
tions in surface gravity turned out to be very important. A ma- 
jor advantage of NSATMOS for this problem was its speed, 
which allowed us to compute extensive grids of models to 
cover the ranges of effective temperature and gravity tailored 
to our needs. Our NSATMOS code will be made available to 
the astronomical community through the XSPEC website 3 . 

Our Chandra observations are described in ^2] We compare 
the models and data in ^3] and discuss the implications in 
|4] Our new neutron star atmosphere models are described in 
detail in Appendix A, and the effects of our consideration of 
varying surface gravity are discussed in Appendix B. 

2. OBSERVATIONS 

The data used in this paper are from the 2000 and 2002 
Chandra observations of the globular cluster 47 Tuc. Both 
sets of obs ervations and t heir i nitial reduction are described 
in detail in lHeinke et alJ (|2Q05|): prior a nalyse s of the 2000 
dataset are described in Grindlav et al. (200lj) and HGL03. 
The 2000 observations were performed with the ACIS-I CCD 
array at the telescope focus, while the 2002 observations 
placed the back-illuminated ACIS-S aimpoint at the focus for 
maximum low-energy sensitivity. Five consecutive observa- 
tions were performed in 2000, as listed in Table 1, with 
three short observations, using a subarray and faster readout 
time, interleaved to reduce pileup in bright sources such as 
X7 and X5 (see HGL03). Pileup occurs when two X-ray pho- 
tons, arriving at the detector during one frame time, are erro- 
neously identified as a single photon with the sum o f the two 
photon energies, or else discarded ("see lDavisll2001l) . Pileup 

3 http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/ 
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TABLE 1 

Summary of Chandra Observations 



Seq, OBSID 


Start Time 


Exposure 


Aimpoint 


Frametime CCDs 


300003,078 


2000 Mar 16 07:18:30 


3875 


ACIS-I 


0.94 


1/4 


300028,953 


2000 Mar 16 08:39:44 


31421 


ACIS-I 


3.24 


6 


300029,954 


2000 Mar 16 18:03:03 


845 


ACIS-I 


0.54 


1/8 


300030,955 


2000 Mar 16 18:33:03 


31354 


ACIS-I 


3.24 


6 


300031,956 


2000 Mar 17 03:56:23 


4656 


ACIS-I 


0.94 


1/4 


400215,2735 


2002 Sep 29 16:59:00 


65237 


ACIS-S 


3.14 


5 


400215,3384 


2002 Sep 30 11:38:22 


5307 


ACIS-S 


0.84 


1/4 


400216,2736 


2002 Sep 30 13:25:32 


65243 


ACIS-S 


3.14 


5 


400216,3385 


2002 Oct 01 08:13:32 


5307 


ACIS-S 


0.84 


1/4 


400217,2737 


2002 Oct 02 18:51:10 


65243 


ACIS-S 


3.14 


5 


400217,3386 


2002 Oct 03 13:38:21 


5545 


ACIS-S 


0.84 


1/4 


400218,2738 


2002 Oct 11 01:42:59 


68771 


ACIS-S 


3.14 


5 


400218,3387 


2002 Oct 11 21:23:12 


5735 


ACIS-S 


0.84 


1/4 


NOTE. — Times in seconds. Subarrays are indicated by fractional numbers of CCDs. 
OBSJD 3385 in this paper, due to its relatively high background. 


We do not use 



has effects upon both spectral and timing analyses, as dis- 
cussed in_the Chandra Proposers' Observatory Guide 4 and 
Davisl (12001 . 

The 2002 observations also interleaved long observations 
with shorter ones to collect some data relatively free of pileup. 
The countrates for X7 in the 2000 and 2002 observations 
(both essentially on-axis) were 0.07 and 0.12 cts/s, leading 
to predicted pileup rates of 9% and 15% for full-array obser- 
vations, or 2% and 4% using 1/4 subarrays. We reprocessed 
the 2000 observations using the CTI correction algorithm im- 
plemented in CIAO 3.2 acis_process_events. Both the 2000 
and 2002 observations were reprocessed (using CIAO 3.2) to 
remove the 0"5 pixel randomization added in standard pro- 
cessing, use updated (time-dependent) gain files, and improve 
hot pixel identifications. Some background flaring occurred 
during parts of the 2002 observations, particularly affecting 
OBS_ID 3385. We do not use data from that short (5 ksec) 
observation, but keep all other data. 

We used the ACIS_EXTRACT software JBroos et alJ 
l200l . version 3.65, to extract and combine spectra and re- 
sponse files for the various observations. We extracted spec- 
tra from contours matching the 95% encircled energy (at 1.5 
keV) of the Chandra point-spread function at X7's position. 
We constructed response files using the mkacisrmf response 
generator in CIAO 3.2. The effective area files take into ac- 
count the decreasing quantum efficiency of the ACIS chips, 
and are corrected to account for the energy-dependent frac- 
tion of the poin t-spread f unction enclosed by the extraction 
region (Broos et al. 2002). We combined spectra from: the 
two long 2000 observations, OBSJDs 953 and 955; the four 
long 2002 observations; the three remaining short 2002 ob- 
servations (see above); and the three short 2000 observations, 
for a total of four summed spectra. We binned these spectra 
at 80 counts/bin for the long 2002 spectrum, 20 counts/bin for 
the short 2000 spectrum, and 40 counts/bin for the other two 
spectra. 



4 http://cxc.harvard.edu/proposer/POG/ 
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2.1. Timing Analysis 

We adjusted all event times to the solar system barycen- 
ter using satellite orbit files provided by the Chandra X-ray 
Center. The qLMXB X5 continues to show eclipses and 
very strong dipping activity throughout the 2002 observations, 
which make a study of its spectrum more complicated; we de- 
fer detailed studies of X5 to a later paper. 

Kolmogorov-Smirnov and Cramer-von Mises tests showed 
no variability from X7 on any timescales probed by the 2002 
observations (seconds to weeks). Power spectra (constructed 
using XRONOS) of X7's lightcurves showed a flat power 
spectru m with less power than expected from Poisson pro- 
cesses (Leahvet al. 1983), typically ~ 70% of the expected 
level. We attribute this to the effects of moderate pileup (as 
other, weaker, X-ray sources in 47 Tuc showed the expected 
levels of white noise), and note that it makes a quantitative 
limit on X7's variability difficult to determine. Assuming con- 
stant Poisson noise (the constant level was a f ree parameter), 
plus re d noise with a fixed slope (cx v~ x , e.g. iRutledge et all 
2001a), we constrain (3a, between 10~ 5 and 0.1 Hz) the rms 
variability to < 3.6%, < 7.3%, < 1.9% and < 5.2% for the 
four 2002 observations. These values may be underestimates 
due to the effects of pileup. We also extract power spectra 
from l-2."annuli, to reduce the effects of pileup, giving 511 
to 541 counts per dataset. These power spectra display white 
noise at the levels expected from Poisson noise. We again 
find no evidence for excess variability, with 3a upper limits 
on rms variability being <17%, <15%, <25% and <28% 
for each observation. The data taken in subarray mode to 
reduce pileup showed similar properties, with 3cr upper lim- 
its of <13%, <17%, and <24% rms variability (excluding 
OBS_ID 3385). No signal was seen at 5.50 hours, the period 
of a marginal signal identified in HGL03, indicating that the 
suggested orbital period (from the lower-quality 2000 data) is 
probably spurious. 

3. SPECTRAL ANALYSIS 
3.1. NSATMOS Spectral Analysis ofX7 

Our standard XSPEC model consists of the NSATMOS hy- 
drogen atmosphere model, absorbed by interstellar gas, con- 
volved with the XSPEC pileup model. We use the XSPEC 
jArnaudll 996) version of the pileup model of J. Davis (Davis 
1200 ll) . setting the frame time parameter to the time resolu- 
tion of each spectrum (header keyword TIMEDEL), and the 
parameters go=l and psffrac=0.95. We floated the grade mi- 
gration parameter a (which parametrizes the fraction of piled 
photons that are recorded as good grades), but found that its 
value is alw ays c onsistent with 0.5, and fixed it to 0.5 for this 
section and 43. 21 

We use the XSPEC phabs model, with lWilms et al.l (120001) 
interstellar element abundances, to describe the interstel- 
lar gas between the Earth and 47 Tuc. We fix its nor- 
malization at Nh = 1.3 x 10 20 c m" 2 , derived using E(B- 
V)=0.024±0.004 as measured by IGratton et al.1 d2003l) . and 
assuming Ry=3.\ JCardelli et al]fl989l) and N H /A v = 1.79 x 
10 21 cm -2 dPredehl & Schmittll 19951) . The Nh column mea- 
sured using any of our models is larger than this, indicating 
additional gas intrinsic to the system (or slight errors in the 
Chandra ACIS calibration below 1 keV). We model this ad- 
ditional absorbing gas by the XSPEC vphabs model, select- 
ing the abundances to correspond to those of 47 Tuc. We 
choose the abundance of iron to be 20% of solar ([Fe/H]=- 
0.7), that of metals between Ne and Ca to be 40% solar 




1 2 5 

channel energy (keV) 

FIG. 1 . — From top, Chandra ACIS-S subarray spectrum (green), ACIS-S 
full-frame spectrum (black), ACIS-I subarray spectrum (blue), and ACIS-I 
full-frame spectrum (red) of X7, fit with our hydrogen-atmosphere NSAT- 
MOS model (histograms) and photoelectric absorption (see text). Pileup is 
responsible for the variation in count rate between the subarray and full-frame 
spectra,and is accounted for in the fit. The apparent difference between the 
full-frame spectra is due to the differing spectral responses of the ACIS-S 
and ACIS-I instruments, as no model parameter is allowed to vary between 
the fits. The strongest residuals (near 2 keV) are instrumental features due to 
the iridium edges of the mirror and resultant rapid changes in effective area 
with energy. See the electronic edition of the Journal for a color version of 
this figure. 



([X/H]=-0.4), C, N, and O to be 63% solar (rX/Hl=-0.2), and 
He at solar abundance (u sing ICarnevl H 996: Salaris & Weissl 
1998: IGratton et al.ll2003l) . We note that our results are not 
very sensitive to the detailed abundances, since the absorption 
column is low (see below). 

The Chandra ACIS effective area and response matrix cali- 
brations are very uncertain below 0.5 keV. We find substantial 
negative residuals around 0.4 keV, similar to an absorption 
line (but substantially narrower than the instrument resolu- 
tion). This feature is s een in most other X-ray sources of suf- 
ficient flux in 47 Tuc (Heinke et al. 2005J, so we ascribe it to 
an instrumental effect and ignore data below 0.5 keV in this 
paper. A second feature is a wave in the residuals between 1.7 
and 2.3 keV. This feature can be ascribed to the difficulty of 
calibrating the iridium M edge structure (Chandra Proposer's 
Observatory Guide), and simil ar waves are seen i n other high- 
quality ACIS CCD spectra (ISanders et alJl2004l) . This resid- 
ual causes the quality of our best fits to be slightly less than 
nominal, Xy=E20 for a null hypothesis probability (nhp) of 
1.5%. More acceptable fits (nhp=7%) can be achieved by 
adding a systematic error term (of order 2%), or by exclud- 
ing the spectral range from 1.9-2.25 keV in the two highest- 
quality spectra, but the parameters of the fit do not change 
substantially from the fits using all the data. 

Fitting the four combined X7 spectra simultaneously in 
XSPEC, leaving no parameters free between the various 
datasets, and only the intrinsic absorption, NS temperature, 
and NS radius as free parameters, gives a reasonably good fit 
(X 2 =l-20, nhp=1.5%; see Fig.[0. (To begin with, we fix the 
NS mass at 1.4 M©. We will vary this later.) The inferred 
(unabsorbed) luminosity of X7 is L x (0.5-10 keV)= 1 .5 x 10 33 
ergs s" 1 , Lb i = 2.4 x 10 33 ergs s" 1 . Allowing the absorption 
column or the NS temperature, radius, normalization, or dis- 
tance to vary between the 2002 and 2000 observations does 
not improve the fit, and the best-fit values for each obser- 
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vation lie within the one-sigma errors of the best fit for the 
other observation. A quantitative constraint on X7's spec- 
tral variability between the 2000 and 2002 observations can 
be obtained by decoupling the 2000 and 2002 temperatures 
or normalizations, and measuring how large the difference 
may be; we find £r X7 ,20oo = l-003^:oo6 x kT xi2002, or for nor- 
malization, K X7 . 2m) = l-01!£o3 x #X7,2002 (90% conf.). Fit- 
ting the individual spectra extracted for each of the four long 
2002 observations (grouped at 20 counts/bin) with the same 
model gives a good fit (xl =0.962, nhp=70%), and when the 
model parameters are allowed to vary between observations 
they again agree with one another. (The improvement in the 
fit quality, compared to the combined spectrum, is due to the 
reduced statistics of each observation individually.) Therefore 
we conclude that these deep observations provide no evidence 
for any spectral change in X7 over a period of 2.5 years. This 
result is in contrast to other qLMXBs which have di splayed 
substantial variability in quiescence (iRufledge et al.l [2002b: 
ICampana et al.l2004l) . 

3.2. Edge or power-law? 

HGL03 found marginal evidence for the existence of an 
edge or other absorption feature near 0.63 and 0.66 keV in 
X7 and X5, respectively. No such feature is apparent in the 
2002 X7 data, and the evidence for such a feature from the 
2000 data has decreased as the Chandra calibration has im- 
proved. Fitting all the X7 data with the standard model above 
and an edge fixed at 0.63 keV (as in HGL03) does not im- 
prove the fit, and leads to a 90% confidence upper limit of 
t < 0.045 for the edge. The two explanations for the edge 
suggested by HGL03 (an OV edge from an ionized wind, or 
a signature of a NS atmosphere that is not purely hydrogen) 
can both be excluded by the disappearance of this feature in 
the 2002 dataset. We feel that the most likely explanation for 
this apparent feature is small uncertainties in the calibration. 

The existence of any spectral feature in a qLMXB spec- 
trum would be of great interest, since an identification of the 
feature could allow determinati on of the gravitational r edshift 
at the NS surface JBrown et alJI 19981 iRutledge et all f2002b. 
HGL03). Therefore we searched for any possible features 
in X7's spectrum, using an edge or a gaussian absorption 
line, between 0.55 and 3 keV. The largest possible features 
were at 0.85 keV, where an edge with r = 0. 065^05 gave 
xl = 1.017; 1.58 keV, where an edge with r = 0.14+oJo gave 
Xl = 1.194; and 2.54 keV, where an edge with r = O.lO^j 
gave xl = 1 -204. The latter two are probably caused by the 
calibration uncertainties around 2 keV. These features have 
F-test probabilities (to be interpreted with caution) of 7.7%, 
13%, and 45% of being generated by chance, which are con- 
sistent with the low significances of their optical depths (none 
nonzero at more than 98% confidence). The 90% confidence 
limits on edge equivalent widths at these locations are < 18, 
< 29, and < 53 eV. Using a gaussian absorption line with 
fixed intrinsic width of 0.05 keV, we find similar results, with 
the most likely features located at 0.92 (r = 0.083![j f 5 ) or 1 .73 
(t = 0.33^jj [5) keV. Equivalent width upper limits are < 1 1 
and < 33 eV, respectively. 

The feature at ~0.4 ke V discussed above is the strongest ap- 
parent feature. The anticipated energy of the strongest feature 
(due principally to the O VIII edge) likely to appear on an ac- 
creting neutron star with a near-solar abundance atmosphere 
is 0.87/(1 +z) keV dBrown et alll998tlRufledge et alJ2 002b. 
see their Fig. 1). There are no plausible features expected near 



0.4 keV for gravitational redshifts in the range 0.2-0.4 (typical 
for favored neutron star equations of state). We conclude that 
no spectral features intrinsic to the NS atmosphere have been 
detected. 

The hard power-law component identified in many 
field q LMXBs may b e a signal of continued accretion 
(Rutledge et al. 2002b), or of nonther mal emission from 
the pulsar wind of an underlying MSP (Bu rderi et alJl2003l 
Bogdanov et al. 2005). We constrain such a power-law com- 
ponent by adding it to our standard model with a fixed photon 
index T of 1.5, 2 or 1, and deriving constraints upon its flux 
in the 0.5-10 keV band relative to the total (absorbed) 0.5-10 
keV flux. Adding this component does not significantly im- 
prove the fit (xl = 1.205, an F-test suggests a 60% chance 
of this level of improvement by chance). For T = 1 .5, such a 
power-law component would make up < 3.2% (90% conf.) of 
X7's total 0.5-10 keV flux. For T = 2 or 1, the constraints are 
< 3.4% of the total 0.5-10 keV flux. 

The complete lack of evidence for variability, edges, or hard 
power-law components in X7 suggests that X7's X-ray emis- 
sion is produced en tirely through re-emission of stored heat 
dBrown et alJll998T) . The constraints on these properties for 
X7, more stringent than for any other qLMXB, make this an 
ideal object for comparisons with NS hydrogen-atmosphere 
models. 

3.3. Constraining the NS Mass and Radius 

The principal goal of our spectral fitting is to self- 
consistently constrain the allowed space in mass and radius 
of X7's neutron star. First we explore the best fits for mass 
and radius if one parameter is fixed at a canonical value, 
M ns =1.4 Mq or /? ns =10 km. Our standard model above keeps 
the mass fixed at 1.4 Mq, and the fitted radius is 14.2+Jq 
km (90% conf.). Allowing the pileup parameter a to vary, 
and allowing for a power-law component with photon index 
fixed at 1.5, the fitted radius is 14.5^j'4 km. The parame- 
ters for this fit are included in Table 2. We note that a 
pure helium atmosphere (possi ble if the donor is degenerate) 
would give even larger radii (IZavlin et al.lll996l iPons et all 
120021). which we th ink unlikely. Adding a 6.1% uncertainty 
JGratton etail2 003. 90% conf., assuming errors are normally 
distributed) in the distance to 47 Tuc gives a NS radius of 
14. 5* j 6 km. The lower limit to X7's radius is thus 12.9 
km, which is signific antly larger than the radius predicted 
by the APR (1 1.5 km.lAkmal et alJI 19981) and FPS (10.8 km, 
Pandharipande & Ravenhall 1989) equations of state for a 1.4 
M Q NS. 

If we instead allow M ns to vary and fix R ns at 10 km, we find 
thatM ns is constrained to 2.17+jj' [^A^; allowing a to vary and 
including a possible powerlaw component, M„ s =2. 20^' [5 Mq, 
Adding distance uncertainty, the errors are M ns =2.20^o Mq, 
Constraining the neutron star mass to lie below the causality 
line gives M ns =2.2O^0^g Mq. This mass is significantl y large r 
than that of any well-measured NS (but cf. iNice et all 2005). 
A high mass for X7 seems improbable given its relatively high 
X-ray luminosity (compared to other qLMXBs), and the ex- 
pected tendency for more m assive NSs to cool more quickly 
( Yakovlev & Pethick 2004). We give the parameters for fits 
with either M ns or R ns held fixed in Table 2, including all un- 
certainties except for the 6% distance uncertainty. 

We use the steppar command in XSPEC to vary both 
the radius and mass parameters, allowing all other free pa- 
rameters (including a and the powerlaw normalization) to 
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TABLE 2 

X7 Spectral Model Parameters 



Model Parameter M fixed R fixed 



Ryhicki NSATMOS modeT 



kT, eV 


105.4+2-5 


151.51" 

4.ot;| 

(10.0) 


Rm, km 


4 2+1 .8 
H,z -1.6 

14.5 + J-« 


M m 


(1.4) 


2 20+°- 16 


xJ/dof 


1.21/251 


1.21/251 


Null hyp. prob. 


1.3% 


1.3% 


Zavlin NSAfc) model 


ifcT.eV 


89+} 


145+ti 


"H20 


5 6+ ' 8 


3 5+ 50 
J - -2.8 


Rns, km 
M„ s 


19.9*}* 
(1.4) 


(10.0) 
1 gg+0.13 




1.21/251 


1.21/251 


Null hyp. prob. 


1.2% 


1.3% 


Gansicke HYD_SPECTRA(#,) model 




z = 0.306 


R fixed 


kT, eV 


123 + | 


138« 


N a 

n H.20 
R U s, km 

M m 
X 2 ,/dof 


c ,+0.8 
J,1 -0.9 
119+1-5 

1 66+°'^ 

1DD -0.15 
1.21/251 


4 4+16 
(10.0) 
1.68*?° 
1.21/251 


Null hyp. prob. 


1.3% 


1.3% 



NOTE. — All en'ors are 90% confidence limits. Distance of 4.85 kpc 
is assumed (its uncertainty is not included in these quoted errors). In each 
column, either mass, redshift (z), or the true radius of the neutron star is 
held fixed. " Nh intrinsic to system, in units of 10 20 cm~ 2 , in addition to 
galactic column of 1.3 X 10 20 cm -2 . * Reached hard limit of model. 




Radius, km 

FIG. 2.— 68% (dotted), 90% (short dash) and 99% (long dash) confidence 
contours in the mass-radius plane derived for X7 by our spectral fitting with 
the NSATMOS model. The causality line, above which no realistic NS equa- 
tions of state can exist, is plotted, along with five representative equations of 
state, APR, FPS, MSI, P CL2, and SQM3 (the last a quark star model, see 
ILattimer & Prakash 2O0jI). 



vary to find the best fit. We show la, 90% confidence, 
and 99% confidence (Ax 2 = 2.3, 4.61, and 9.21 respec- 
tively) contours in NS mass and radius in Fig. [2] for X7, 
using our spectral model as described above. These can 
be compared with the l oci of models following the APR 
JAkmal et alJ 19981). FPS llPand harinand e & Ravenhal]ll989l) . 
PCL2 dPrakash et alJl995])7 MSl (Miiller & Sero tl 19961) . and 
SQM3 ( Prakash et all 19951) equations of state for dense mat- 
ter (the last is a representative "quark star" model). For the 
SQM3 and APR models, a high NS mass is required for con- 
sistency with the data; the FPS model is ruled out; and the 
PCL2 model is marginally consistent with the data at the 99% 
level for a mass of 1.35 M . The MSI model is an exam- 
ple of models consistent w ith a larger radius at 1 .4 M Q (see 
Lattimer & Prakash 2001), some of which include condensa- 
tion of kaons, hyperons, muons or free quarks in the core. 
The line labeled "Causality" (R = 3.04GM/c 2 ) represents the 
requirement that the speed of sound must be less than c, a 
necessary (but not sufficient) requirement for a NS interior 
to respect causality (Lattimer & Prakash 2001; Olson 2001). 
This requirement is virtually identical to R > 3GM/c 2 , which 
is the condition that the neutron star surface is outside its own 
photon sphere, so that self-irradiation does not occur. 

3.4. Comparison ofNS atmosphere models 

We compare three different hydrogen-atmosphere neutron 
star models in this analysis. Our NSATMOS model has pa- 
rameters M ns (mass of neutron star, M ), R ns (true radius of 
the neutron star, km), T e g (the effective temperature of the 
neutron star surface, unredshifted), and D (distance to neutron 
star, pc). The NSA model (Zavlin et al. 1996) has parameters 
R (true radius of the neutron star, km), M (mass of neutron 
star, M Q ), Teff (the effective temperature of the neutron star 
surface, unredshifted), and normalization K = 1 /D 2 (where D 
is the distance to the neutron star, pc). The last parameter is 
often used to calculate R^ = R(l + z), the radius as seen by 
a distant observer, since the last parameter can be considered 
K = (Roo/Roo,o) 2 /D 2 , where Roo.o is the equivalent R^ that 
would be calculated from the mode l parameters R and M. The 
HYD_SPECTRA model of Gansic keetail J2002T) has the pa- 
rameters r e ff, z (the gravitational redshift at the NS surface), 
and normalization K = (R/D) 2 , where R is in units of 10 km 
a nd D is measu red in pc. We also tested the H_ATM model 
of Llovd (2003), for which the parameters are r e ff, z, and nor- 
malization K = (R/D) 2 , where R is in km and D is in units of 
10 kpc. The H_ATM model gives results broadly similar to 
those of the NSA model, and for brevity we will not discuss 
H_ATM further in this paper. The NSATMOS and NSA mod- 
els fit directly for M ns and R NS , while the HYD_SPECTRA 
model fits redshift and normalization, which allows computa- 
tion of M ns an d ^ns- 

While the NSA and H YD_S PECTRA codes can, in princi- 
ple, be used to compute models for any value of gravity, the 
actual versions of these models found in the standard XSPEC 
package have been constructed for single values of the sur- 
face gravity g s (appropriate for a 1.4 M Q , 10 km NS). To em- 
phasize this, we shall denote these models by NSA(gj) and 
HYD_SPECTRA(g J ). The use of these fixed gravity mod- 
els for other values of surface gravity is not strictly appro- 
priate. In order to account for the effect of surface gravity 
on the spectra, we computed an extensive grid of models us- 
ing the NSATMOS code described in Appendix A. This grid 
was used as the basis of an interpolation routine that produced 
model spectra for the XSPEC program for a full range of ef- 
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FIG. 3. — 90% confidence contours in the mass-radius plane derived for X7 
by our spectral fitting with hydrogen atmosphere NS models: N SATMOS, us- 
ing variable g, with thick solid (blue) lines, NSA(g t ) i Zavlin et al. 1996) with 
dashed (green) lines, HYD_SPECTRA(g,,) I Gansicke et alJ2002l) withshort- 
dashed (red) lines, and, for comparison, NSATMOS(g J ), with thick dotted 
(magenta) lines. Here g s = 14.385 is the surface gravity for the standard NS 
model. See the electronic edition of the Journal for a color version of this 
figure. 



fective temperatures and gravities. 

We fit each of these three models to our X7 data using the 
same spectral model as described in section 33.31 replacing 
only the neutron-star atmosphere model. For constraints upon 
the NS mass, radius, or gravitational redshift, we place the 
parameters of these fits in Table 2. We plot the contours en- 
closing 90% confidence contours for various models together 
in Fig. |3] The contours labelled NSATMOS were constructed 
using the appropriate surface gravity for each point in the 
plot. We have also plotted the contours assuming a fixed stan- 
dard gravity logg., = 14.385 for three codes: NSATMOS(g. 5 ), 
NSA(g s ), and H YD_S PECTR A(g s ) , the latter two using the 
models in the XSPEC package. 

We see that none of the available models allow a 1 .4 M , 10 
km radius NS, requiring a more massive or larger radius NS. 
[However, the small additional uncertainty on the distance 
to 47 Tuc-not included in Fig. [3]-allows marginal consis- 
tency at the 99% level with the HYD_SPECTRA modelfe).] 
Since the HYD_SPECTRA(g. s ) and NSA(g. s ) models are con- 
structed with a g s appropriate for a 1.4 M Q , 10 km NS, the 
failure of those models to find an acceptable fit for a 1.4 
Mq, 10 km NS is robust. It is not clear to us why the 
HYD_SPECTRA(gj) model predictions are significantly dif- 
ferent from the other two models, although the difference be- 
tween the HYD_SPE CTRA and NSA hydr ogen models has 
been previously noted (Gansicke et al. 2002). Caution should 
still be taken in interpreting our results due to possible un- 
known systematic uncertainties in the aAChandra calibration. 
However, the robustness of all tested models in excluding the 
canonical mass and radius values ( 1 .4 M Q and 1 km) is strong 
evidence for a relatively "stiff" NS equation of state. 
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FIG. 4. — 90% confidence contours in the mass-radius plane derived for 
X7 by our spectral fitting with two versions of our hydrogen atmosphere NS 
models: NSATMOS (variable g) and NSATMOSfc), compared to lines of 
constant Roo. See the electronic edition of the Journal for a color version of 
this figure. 



4. DISCUSSION 

In Fig. |4j we compare the allowed contours of mass and 
radius for X7 using the full NSATMOS and fixed-gravity 
NSATMOS(gs) models to lines of constant R^ for a black- 
body spectrum (or for any spectrum of the scaled form in 
equation [Bl]). The contours derived from our grid of NSAT- 
MOS models lie closer to lines of constant R x than do the 
contours derived from the fixed-gravity models, but there are 
significant differences. (See Appendix B for further discus- 
sion of the surface gravity effects.) 

The effect of self-consistently including variation in g s is 
quite substantial, as can be seen from the difference between 
our NSATMOS contours and the contours of NSATMOSfc). 
The predictions of the NSA(g s ) model are very similar to 
those of our NSATMOS(g s ) model, demonstrating that vari- 
ation in surface gravity is the primary reason for the dif- 
ferences in the results obtained with NSATMOS and NSA. 
This demonstrates that the NSA(g. s ) and HYD_SPECTRA(g. s ) 
models currently available in XSPEC do not give correct con- 
straints in regions of parameter space where the assumed sur- 
face gravities are different than those for which they were 
computed, and which are outside their region of validity. Ex- 
plicitly, the constraints on neutron star mass and radius com- 
puted by HGL03 used these models outside their range of va- 
lidity, and produced incorrect results. It is possible that this er- 
ror in assumptions may affect the results of lPons et alJ d2002h 
and lWalter & Lattimerj j2002) as well. 

Some works have us ed the NSA models in XSPEC to con- 
strain th e range of (Rutle dge et alJ2001blaHGendre et alJ 
2003©. This method is more accurate due to a degener- 
acy in spectral shape vari ations between surf ace gravity and 
surface temperature (see Zavlin et al. 1998, and Appendix 
B). For instance, fitting X7 with the NSA model, freezing 
the mass (1.4 M Q ) and radius (10 km), we infer R^ = 18.34 
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km, with a 90% confidence range between 17.1 and 22.1 km. 
Comparing this to Figure 4, we see that this range is more 
accurate than the method of HGL03, though still inaccurate 
at the ~10% level. The differences between hydrogen atmo- 
sphere spectra and blackbodies are large enough that direct 
calculation of a range of M ns and R ns values giving acceptable 
fits is preferable to simply constraining Roo . Measuring R^, is 
generally satisfactory, however, when the observational data 
is not of the highest quality, and/or the purpose of the mea- 
surement is simply to check for consistency with the canoni- 
cal NS mass and radius (a typical test of current observational 
data). 

Applicability of any pure hydrogen-atmosphere model to a 
NS atmosphere requires that accretion is not continuing above 
the critical rate to keep metals present in the atmosphere, as 
that w ould alter the opacity of the atmosphere iBrown et alJ 
QUI). TheX -ray luminosity produced by the critical accre- 
tion rate is similar to the observed X-ray luminosity of X7. 
Thus, it remains conceivable that the inferred mass and radius 
of X7 are biased by the presence of metals, and thus extra 
opacity, in the NS atmosphere. However, the extraordinary 
stability of X7's X-ray flux, and its lack of accretion signa- 
tures such as a hard power-law spectral component, or fea- 
tures due to lines or edges, indicate that X7's X-ray emission 
is most likely produced by deep crustal heating rather than 
continued accretion. 

5. CONCLUSIONS 

We have computed new grids of hydrogen atmosphere mod- 
els for neutron stars in quiescent LMXBs, accounting self- 
consistently for electron thermal conduction, radiation force, 
self-irradiation of the NS, and variations in the NS surface 
gravity. The first three effects do not produce large changes, 
and our NSATMOS models agree very well with the NSA 
models of ( Zav lin et all 199 6) in the appropriate surface grav- 
ity regime. We find that the effects of freely varying the NS 
surface gravity, on the other hand, are significant, especially 
when trying to constrain the mass and radius of neutron stars. 

We report new (2002) Chandra observations of the quies- 
cent LMXB known as X7 in 47 Tuc. No variability is ob- 
served on any time scale between minutes and weeks in the 



new data set. The 2002 Chandra data, plus the 2000 Chandra 
data, can be simultaneously fit with our hydrogen atmosphere 
model, photoelectric absorption, and a correction model for 
instrumental pileup. No convincing spectral features are seen, 
and no additional hard component is detected. No variations 
are seen in X7's spectrum over the 2.5 year interval. 

Our spectral fitting constrains the range of mass and radius 
which can produ ce a spectrum like that observed. In contrast 
to the results of Heinke et al. (2003), our use of a range of 
surface gravities allows acceptable fits to X7 with a standard 
1.4M Q NS mass and a radius of 14.5*}" ® km, as well as a high- 
mass solution (M=2.20!J]; ( ^ M ) for a radius of 10 km. For 
a canonical mass of 1.4 Mq a 10-12 km radius is ruled out 
at 99% confidence. This indicates (assuming the validity of 
a pure hydrogen atmosphere model) that either X7 is more 
massive than any NS yet known, or that X7 has a somewhat 
larger radius than canonical modern NS models (our preferred 
interpretation). In either case a relatively "stiff" NS equation 
of state is favored. 

The HYD_SPECTRA ( Gansi ckTeTafl l200l and NSA 
hydrogen-atmosphere NS models as currently implemented 
in XSPEC do not include a range of surface gravities appro- 
priate for the possible ranges of NS masses and radii. Our 
work using NSATMOS has shown that accurate constraints 
(at the <10% level) on the radius and mass of NSs require 
self-consistent modeling of the effects of surface gravity. Our 
NSAT MOS model, and the NSAGRAV code (using the mod- 
els of lZavlin et all 19 96. provided for XSPEC during this pa- 
per's refereeing process), meet these requirements. 



We warmly acknowledge Don Lloyd for providing compar- 
ison neutron star atmosphere models for use in verifying our 
code, and Marc Freitag for the use of his fig2curve script. We 
also thank Lars Bildsten, Bob Rutledge and Jim Lattimer for 
useful discussions and suggestions, and the CXC team for 
their untiring data calibration efforts. This work was sup- 
ported in part by NSF grant AST 0307433, and in part by the 
Lindheimer Postdoctoral Fellowship at Northwestern Univer- 
sity. 



APPENDIX 
THE NSATMOS CODE 

The NSATMOS code calculates an atmospheric model under the following conditions and assumptions: (1) Static, atmosphere 
in the plane-parallel approximation; (2) Negligible magnetic fields {B < 10 8 G); (3) Ideal equation of state for pure hydrogen 
with complete ionization; (4) Opacity due to thermal free-free absorption plus Thomson scattering in the unpolarized, isotropic 
approximation; (5) Energy transport by radiation and electron heat conduction; (6) Hydrostatic equilibrium includes radiation 
force due to absorption and scattering; (7) Comptonization is ignored; and (8) For a compact object within its own photon 
sphere,/? < (3/2)Rs, self-irradiation of the surface is taken into account (Rs, = 2GM/c 2 is the Schwarzschild radius). We note 
that the neglect of neutral hydrogen limits the validity of the code to te mperatures T e ff > 3 x 10 5 K, while the omission of 
Comptonization limits it to perhaps T e ff < 3 x 10 6 K dZavlin et alJl!996l) . These limitations are not serious for the present 
application, but work is in progress to correct NSATMOS for both. 

The solution of the atmosphere problem based on these assumptions follows mostly along standard lines for stellar atmospheres. 
Here we review the main ideas for reference. 

The depth variable used here is the mass column density m, measured from the surface, which is related to vertical height z by 
the differential relation dm = -pdz, where p is the mass density. 

The intensity field 7„(m, p) is a function of frequency v, depth, and p, the cosine of the angle relative to the outward normal. The 
transfer equation for the intensity is then 

dl 

p-^ = K v {I v -B v ) + K T (I v -J v ). (Al) 

The Planck function is denoted by B = B U {T), where T = T(m) is the local temperature. For a fully ionized hydrogen plasma the 
Thomson opacity kj is a constant, but the free-free opacity n u depends on the density and temperature. The mean intensity J u is 
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defined by the integral, 

i r i 

H J Udfi, (A2) 



2 



for j = 0. Other useful moments H u and K v are defined by this integral for j = 1 and j = 2, respectively. 

When the stellar surface is outside the photon sphere, R > (3 /2)Rs, the surface boundary condition on the intensity field 7„(m, /j) 
is the usual one of no incident radiation. When R < (3 /2)R$, the self-irradiation of the surface due to gravitational bending of the 
rays is taken into account as discussed in McClintock et al. (2004). The boundary condition at suffiently large depth is that the 
radiation field is given by the LTE diffusion approximation based on the local gradient of the temperature field. 

An important quantity is the radiative monochromatic energy flux, given in terms of the H„ moment by F v = AttH u . The total 
radiative energy flux is found by integration over frequency, 

/>oo />oo 

F mi = / F v cLu = Att I H v dv. (A3) 
Jo Jo 

For some ranges of parameters of interest, the heat conduction by electrons can be of importance. For the conductive flux we use 
the Spitzer-Harm formula 

^on = -A c r 5 / 2 ^, (A4) 
oz 

The constant A c has a standard value of order 1.8 x 10~ 5 erg cm" 1 s" 1 K -7 / 2 . Magnetic field s tend to suppress the thermal 
conduction, though by only a factor of order 3-5 when the field is chaotically tangled (Naravan & Medvedev 2001). Therefore, 
for the calculations of this paper we have assumed a conductivity equal to a third of the Spitzer value. However, we find that 
electron conduction even at the Spitzer value is a negligible contributor to the energy flux, making the issue moot. We take as a 
boundary condition at the surface that the conductive flux vanishes, F con (Q) = 0, since no electrons exist above the outer surface. 
Energy equilibrium requires that the net outward energy flux, radiative plus conductive, be constant with depth, 

F t ot = F m d + F con = a r e ff 4 . (A5) 

The constant total flux is parametrized by the effective temperature T e ff. Since the conductive flux vanishes at the surface, the 
effective te mper ature is, as usual, a measure of the total emitted radiative flux. 

Equation (IA5i is used in NSATMOS as the basic expression of energy equilibrium. We note that another common method of 
expressing this is through the flux derivative condition, 

^=0. (A6) 
dm 

If this form is used, the desired value of total flux F tot must be introduced in some other way, usually through the boundary 
conditions at depth on radiation and conduction. 
The hydrostatic equilibrium equation is 

dp 

Q^=g-g™d, (A7) 

where p is the gas pressure and g is the local acceleration of gravity, given by the usual formulas in Schwarzschild geometry. 
The radiative weakening of gravity g m d is given by the radiation force per unit mass due to free-free absorption and Thomson 
scattering, 

gmd = - / (k» + Kt)F v dv. (A8) 



c t 

In a completely ionized hydrogen plasma, the electron and proton densities are equal, n e = n p , and the ideal equation of state is 

p = 2n e kT (A9) 

where k is the Boltzmann constant. 

The model atmosphere problem requires the simultaneous, self-consistent solution of the preceding equations for the radiation 
field and the gas properties (temperature, density, and pressure) as functions of depth. Overall, this is a set of nonlinear, coupled 
equations. Our method of solution involves (as do all such methods) a preliminary recasting of the equations into favorable forms 
such that cycles of iteration between them converge reasonably rapidly to the desired solution. 

One very useful trick is to recast the transfer equation in terms of Eddington factors. Multiplication of equation (IA1> by 1 and 
fi, followed by integration over all /i gives 

—^ = e v {J v -B v ), (A10) 

^=H„, (All) 
where the monochromatic optical depth is defined differentially as <f r„ = (k v + kj) dm with t„ = at the surface, and where 

6= Kv . (A12) 
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Introducing the Eddington factor f v = K v jJ v , equations JA10> and IA1H can be expressed as the single second-order equation, 

d 2 (fM 



drl 



■ e v {J v -B v ), (A13) 



The two boundary conditions on this equation at the surface and at depth are expressed through the boundary Eddington factors 
guO and g ud , such that, 

H v {0) = g v0 J v {0), H v {T vd ) = g vd J v {T v d). (A14) 

After a set of initial Eddington fac tors is chosen, they are updated during the course of iterative solution using full angular formal 
solutions of the transfer equation \A\\ . It is during this formal solution that the detailed boundary conditions on intensity come 
into play, including the possibility of sef-irradiation. 

A particularly critical part of the iterative solution is how the run of temperature with depth T(m) is corrected after each iteration. 
In doing this it is important to take account of how radiation can couple the variables at distant points in the atmosphere. Methods 
that try to correct the temperature based on local information, or even on information just a few optical depths away, will converge 
much too slowly to be practical. 

Here we adopt a temperature correction scheme based on a partial linearization of the temperature dependence in the transfer 
equation and the energy equation. The temperature field is represented as 

T(m) = T (m (m) + T (l) (m), (A15) 

where T (0 \m) is an initial "guess" for the temperature, and T (l \m) is the "correction." In the perturbation sense these are regarded 
as of zeroth and first order, respectively. Other quantities in the problem can be expanded to first order in the perturbation. For 
example, the radiation field has the representation, 

Wf+e (A16) 

Not all quantities in the problem are expanded in this way, making this a partial linearization rather than a complete linearization 
scheme. The decision as to which quantities are expanded and which are not is somewhat arbitrary, but is guided by the physical 
idea that the critical equations are the energy equation, the transfer equation, and the conduction equation, since these are the 
ones that control the global redistribution of ener gy in the atmosphere. 
The zeroth order form of the transfer equation (IA13b is 

d -^^l = e v [jf ) -B v (T^)]. (A17) 

We take as the first order form, 

d 2 (fvJ ( J } ) _ _ r/ (l)_n / r (0K T (l) 



drl 



eAJi l, -BAT m )T w ], (A18) 



where B(T) = dB u (T)/dT . Note that the radiation fields and the temperature dependence of the Planck function have been 
expanded, but not the Eddington fact ors o r the opacities, so that the optical depth scale is unchanged. 
Expansion of the energy equation iA5\ through first order gives, after some rearrangement, 

*8? = F%+F£l = ^ eS 4 - (JS = -E. (A19) 

The right hand side is the negative of the flux "error" E in the zeroth order solution. The left hand side is the total first order flux, 
consisting of the two terms, 

/•OC 

*S = 47r / H^dv, (A20) 



Jo 

1 con "C q 



(•j(0)n5/2j(1) 



(A21) 



Each of these depends linearly on the first order temperature correction T (r K For F^l this follows directly from equation (IA2H . 

H$\ (A22) 



For ,Fj ( ad we note that depends linearly on since 



(i) . 

wc iiulc LiiciL ii ucjjciiua uncoil^ un , ainv^c 



3t v 

which follows from equation JA1 1> . noting = f v J^. Finally, fjp depends linearly on through equation (IA18I 
In terms of Nd discretized depths, the linear dependence of F^(m) and F^(m) on r (1 ^(m), can be expressed as 

F^ = ^-T m , (A23) 
F^ = <f, con .T m , (A24) 
where $ ra d and $ con are finite matrix operators over the discrete depths. 
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In orde r to c onstruct the matrix $ ra d, Nf discrete frequencies are introduced. For each of these frequencies, the transfer 
equation 1A17I . plus its boundary conditions, becomes a tridiagonal matrix equation in depth, which may be solved in of order 
Nf) operations. In this way one obtains matrices A v such that j£ ' = A v ■ for each frequency. Implementing equation (IA22> with 
some form of numerical differentiation, we obtain matrices B v such that Fjp = B v ■ T^ l \ Implementing a frequency quadrature 
scheme for equation JA20> . we then obtain the matrix $ ra( j. The whole process takes of order NfN^ operations. 

The co nstruc tion of the matrix $ con is easier, since it merely involves implementing a numerical differentiation formula for 
equation JA2H plus boundary condition. Depending on the order of the formula, the resulting matrix might be bidiagonal or 
tridiagonal. 

Defining $ tot = <i> ra d + <i> con , equations ( IA19l l. JA23L and ( IA24> can be combined to give, 



$ tot -r (1) =-£, (A25) 

which relates the temperature correction directly to the flux error E in the trial (zeroth order) solution. This matrix equation 
is then solved for and a new trial temperature law is obtained from T (0) <— T (0) + 

Using the new trial temperature, improved radiation fields, Eddington factors, radiative force, and hydrostatic equilibrium are 
computed. This process is repeated until convergence is obtained. 

This method was used to compute our grid of models. In order to initialize the first model, we used a Rosseland-type solution 
for the run of temperature (similar to the asymptotic result given in equation (15) of Zavlin et al. 1996) and we set the radiative 
force to zero. For subsequent models we adopted those quantities from a previously calculated model with nearby parameters. 
This was particularly helpful in computing models close to the Eddington limit. 

Some additional "tricks" we found useful for some models were: (1) Adiabatic turning on of parameters, in which the model 
parameters were changed in small steps from those of a previously solved case to the desi red on es, these changes occuring along 
with the other iterations; and (2) Linearization limiters, which are virtually identical to (IA15I for small corrections, but which 
keep the temperature correction within fixed bounds, for example, 



T{m) = J (0) (m) + r< tanh 



2T m (m) 



T<®(m) 

SURFACE GRAVITY EFFECTS 



(A26) 



Our spectral fitting to constrain X7 in the mass-radius plane has shown the importance of considering neutron star models with 
a full range of surface gravities as well as effective temperatures. This conclusion is based on numerical results of the XSPEC 
analysis program, which involves a long chain of data reductions plus manipulations of the the model atmosphere spectra. In this 
appendix we attempt to understand better the importance of the surface gravity g s using some heuristic, qualitative arguments. 

We begin by considering a approximate analytical representation for the surface spectral flux of the NS atmosphere, namely, 

F„(v) = T eB 3 <f>(Q, Z = v/T eS , (Bl) 

where (f>(£) is a fixed function. This implies that there is a single spectral shape which is simply rescaled by the effective 
temperature r e ff. Such a representation would hold exactly if the neutron star surface radiated as a blackbody. It also holds 
approximately for more realistic pure ionized hydrogen atmospheres; an approximate expression of this form was given in 
equation (A 17) of McClintock et al. (2004). 

Given the spectral flux at the surface, the observed spectral flux FobsC^obs) at a large dis tanc e D may be expressed (see, e.g., 
McClintock et al. 2004). Let us assume for the moment that the surface flux is given by IB It exactly. Then for a neutron star 
with mass M m and radius R m we may write, 

R 2 

F bs{v ohs )= ^-T o'(j)(v ohs /T 00 ), (B2) 

where = T e g/(1 +z), R^ =R ns (l +z). Here 1 +z = (1-Rs/RnsT 1 ^ 2 is the gravitational redshift factor at the neutron star surface, 
and R& = 2GM as /c 2 is the Schwarzschild radius. 

With a prescribed form of the function <fi and a known distance D, spectral fitting is reduced to the determination of just two 
parameters, Too and R^. However, there are three parameters required to specify a neutron star model fully, which can be taken 
to be its mass M ns , radius R m , and effective temperature T. This implies that for each fit there is a one-parameter set of models 
that are observationally indistinguishable from that fit. To show this in detail, we first note that the relation R^ = R ns (l +z) can 
be expressed, 

M ns= 1 7 -T^T- ( B3 ) 



Roo 2 J 2G 

An acceptable M ns -7? ns pair can lie anywhere on one of the contour lines of constant R^, shown as the solid curves in figure 
IBlT a). The additional parameter along each curve can be taken to be the ratio T^/Too, which is equal to the redshift factor (1 +z), 
given in terms of R ns and M ns by 

2GM n x " 1/2 



1«-U-3==I • (B4) 



Curves of constant (1 +z) = T^/Too are plotted as dashed curves in figure lBTT a). However, since there is no independent determi- 
nation of r e ff, this relation does not constrain the possible locations along the contours of R^. In a spectral fitting program such 
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(a) 3 
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(b) 3 
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FIG. B 1 . — (a) The M ns -i? ns plane, showing solid curves of constant and dashed curves of constant (1 +z). Two of these dashed curves are labelled for 
the values log(l +z) = 0.04 and 0.14. The others are at logarithmic spaced values with Alog(l +z) = 0.02. (b) The same, except with dashed curves of constant 
surface gravity g s . The two marked curves are for logg s = 13.6 and 14.4, with other curves spaced with Alogg s = 0.1. 
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FIG. B2. — (a) Scaled surface number flux vs. scaled frequency for logg s = 14.4 and for logoff ranging from 5.9 to 6.4. (b) Scaled surface number flux vs. 
scaled frequency for logoff = 6.2 and for logg s ranging from 13.8 to 14.8. 



as XSPEC, this degeneracy would sh ow i tself by the fact that the confidence contours would be made up of contours of constant 
Roo, that is, the solid curves in figure lBlT a). 

In practice, since the simple scaling law JB 1 i is only approximately valid, the confidence contours will not follow those solid 
curves precisely. Nonetheless, one can see in figures |5] [3] and0]that the confidence contours are generally consistent with there 
being a blackbody-like degeneracy along their length, since they do not close off on either ends of the confidence bands. One 
also sees that the shapes of the confidence contours are different depending on whether the gravity has been fixed or allowed to 
vary freely. The fixed gravity curves can be characterized as being pushed out to larger radii at smaller mass, giving the curves in 
this region a flatter appearance. The curves for freely varying gravity do the opposite: they are pushed to smaller radii at smaller 
mass, so much so that they actually curl around and change their direction. (It is this fact that allows a fit with a smaller radius, 
one of the major results of this paper.) 

We can give a heuristic, qualitative argument for this difference. First, let us note that the surface gravity, which is required for 
the construction of the atmosphere, is given by, 

gs = —a ■ (B5) 

In figure lBlT b). curves of constant (solid) are again plotted in the M ns -/? ns plane, but now wi th cu rves of constant g s (dashed). 

The sensitivity of the function <f>(v/T e ff) to temperature and gravity is demonstrated in figure lB2l Here is plotted the quantity 
log(/v/i/r e ff 2 ) versus log(^/r e ff) for the surface fluxes of the neutron star models for various effective temperatures lB2f a)1 and 
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various gravities (B2jb)] in the neighborhood of the approximate fitted values log T e s = 6.2 and logg s = 14. 4. Th e frequency range 
here was chosen to match roughly the portion of the flux curves observed in X7. If the scaling relation JB lb were exact, these 
would all lie on the same curve, namely, log[(/>(£)/£]. In fact, there is a sensitivity to both effective temperature and gravity. We 
see that the flux distribution steepens for increasing temperature, and for decreasing surface gravity. The sensitivity to temperature 
is greater than that of gravity, in the sense that 0.5 dex variation of tempe ratur e give s a s ubstantially greater variation than 1.0 
dex variation in gravity. However, an essential point to notice from figures lBTT a) and lBlf b) is that over the part of the M ns -R ns 
plane where the spectral fitting results are most divergent, say, for 0.5 < M ns /M < 2 and for 12 km < R m < 18 km, the gravity 
changes by a bout 1 .0 dex, while the temperature changes by about only 0. 1 dex [a range corresponding to two neighboring curves 
in figure 123 a)]. In that case, the change due to gravity is effectively more substantial than that due to temperature. Also, the 
changes in the curves are of the opposite sense as one moves along the curves of constant R^, from the top to bottom of the plot, 
since both temperature and gravity decrease. Thus, if one fixes the gravity, the confidence contours are pushed in one way by the 
sensitivity to temperature, but if gravity is allowed to vary, then the larger effective sensitivity to gravity pushes the confidence 
curves in the opposite direction. The preceding heuristic argument accounts qualitatively for the overall behavior seen in figures 
|2j[3] and|4]and demonstates the need to include variable gravity when doing spectral fitting for such objects. 
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